Clustering by mixing flows 
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We calculate the Lyapunov exponents for particles suspended in a random three-dimensional flow, 
concentrating on the limit where the viscous damping rate is small compared to the inverse corre- 
lation time. In this limit Lyapunov exponents are obtained as a power series in e, a dimensionless 
measure of the particle inertia. Although the perturbation generates an asymptotic series, we ob- 
tain accurate results from a Pade-Borel summation. Our results prove that particles suspended in 
an incompressible random mixing flow can show pronounced clustering when the Stokes number is 
large and we characterise two distinct clustering effects which occur in that limit. 
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This letter describes the dynamics of particles sus- 
pended in a randomly moving incompressible fluid which 
we assume to be mixing: any given particle uniformly 
samples configuration space. At first sight, it seems as if 
the particles suspended in an incompressible mixing flow 
should become evenly distributed. This indeed happens 
if the particles are simply advected by the fluid. How- 
ever, it has been noted |l| that when the finite inertia 
of the suspended particles is significant, the particles can 
show a tendency to cluster. 

The current understanding of this remarkable phe- 
nomenon refers to a dimensionless parameter termed the 
Stokes number, St = 1/(77"), where 7 is the rate at which 
the particle velocity is damped relative to that of the fluid 
due to viscous drag, and r is the correlation time of the 
velocity of the fluid. There is a consensus 0, IE El IE @ 
that clustering is most pronounced when St is of order 
unity. 

In this letter we argue that strong clustering can occur 
when St is large. We show that different clustering mech- 
anisms compete at large values of St and quantify under 
which circumstances clustering occurs. Before describing 
our results and outlining how they are derived, we briefly 
summarise previous theoretical work on the clustering of 
inertial particles in turbulent flows. 

This effect was first discussed by Maxey Q: he ap- 
proximated the inertial particle dynamics by advection 
in a 'synthetic' velocity field which was obtained as a 
perturbation of the velocity field of the fluid, u(r,t). 
Maxey showed that this synthetic velocity field has neg- 
ative divergence when the vorticity of u(r,t) is high or 
its strain-rate low, and predicted that particles would 
have low concentrations in regions of high vorticity due to 
this 'centrifuge effect'. This effect has been demonstrated 
in direct numerical simulation of particles suspended in 
a fully-developed turbulent flow [E • The theoretical 
work of Maxey and experimental work on turbulent flows 
|E| has emphasised instantaneous correlations between 
vortices and particle-density fluctuations. 

Later work has adapted results on the density statis- 



tics and Lyapunov exponents of purely advective flows 
obtained in Elperin 0] suggested combining these 

results with Maxey's synthetic velocity field to obtain re- 
sults for inertial particles; a similar approach was used in 
0,0, El- These results are not applicable at large St, 
because the perturbation of the velocity field need not be 
small when inertial effects are important. 

An alternative viewpoint arises from work of Sommerer 
and Ott 0], who describe patterns formed by particles 
floating on a randomly moving fluid. They characterise 
the patterns in terms of their fractal dimension and sug- 
gest that the fractal dimension can be obtained from ra- 
tios of Lyapunov exponents of the particle tra ject ories 
using a formula proposed by Kaplan and Yorke 0] . 

The argument in [14j extends to particles suspended in 
turbulent three-dimensional incompressible flows. Con- 
sider the Lyapunov exponents Ai > A2 > A3. They are 
rate constants defined in terms of the time dependence 
of, respectively, the length 5r of a small separation be- 
tween two trajectories, the area 6 A of a parallelogram 
spanned by two separation vectors and the volume 5V of 
a parallelepiped spanned by a triad of separations: 



Ai = lim t 1 \og c (5r) 

t — >OG 

Ai + A 2 = lim rHog^SA) 

t — »oo 



Ai + A 2 + A3 = lim t 1 logJSV) 

t — >oo 



(1) 



The Kaplan- Yorke estimate for the fractal dimension in 
a three-dimensional incompressible flow is determined by 
the dimensionless quantity ('dimension deficit') 



A = -(A 1 +A 2 + A 3 )/|A 3 | 



(2) 



When A > 0, the Kaplan- Yorke estimate of the dimen- 
sion is dn = 3 — A, and dn = 3 if A < 0. Clustering 
effects are significant if the fractal dimension is signifi- 
cantly lower than the dimension of space. This proposi- 
tion provides a strong motivation to study the Lyapunov 
exponents of the problem. 
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A third mechanism for clustering is the following: 
nothing prevents the infinitesimal volume element SV 
from collapsing to zero for an instant of time. These 
events correspond to 'caustics', where faster moving par- 
ticles overtake slower ones. Caustics are associated with 
the density of particles on a surface becoming very high, 
facilitating the aggregation of suspended particles. This 
mechanism was recentl y p roposed as a cause of cluster- 
ing of inertial particles |lfij | , and is also mentioned briefly 
in 0] ■ The significance of this effect is determined by 
the rate J at which the infinitesimal volume element goes 
through zero for a given triplet of nearby trajectories. 

Which of these three mechanisms is most important? 
Maxey's centrifuge effect is weak at small St, where the 
particles are simply advected. There is a consensus that 
the effect is also weak for large St, because the vortices 
do not persist for a sufficiently long time to be effective, 
implying that significant clustering is only observed when 
St ~ 1. However, there is at present no understanding of 
what happens at large values of St. In the following we 
describe quantitative results for the Lyapunov exponents 
Xj , for the dimension deficit A and for the rate of caustic 
formation J: these are summarised in Fig. 1 a - c. 

Our results show that in order to understand the clus- 
tering effect it is necessary to consider not only the Stokes 
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FIG. 1: a Lyapunov exponents as a function of e ~ kxj -1 / 2 : 
results of simulations of an Euler discretisation of the lin- 
earised equations of motion @, replacing the elements of 
F(t) by random numbers with appropriate statistics, sym- 
bols. Also shown are results based on Pade-Borel summation 
of 1171 . with Pade approximants of the orders P22 (— • — ), 

P23 (— ■ — ), P33 ( ), and P34 (solid lines), b shows A vs. 

e (obtained from using the data in a), c shows — -W7 
(circles, solid red line, see a and b) compared to the rate of 
caustic formation determined from simulations as described 
above (squares) and a fit to C exp(— S/e 2 ) with C = 0.21 and 
S — 1/8. d schematic phase diagram in the re - uj plane: 
note the logarithmic scales. The red line indicates where the 
quantity A is zero. Above the red line A is positive implying 
clustering. In the under- and overdamped advective limits A 
positive and small (weak clustering), while below the red line 
A is negative. In this regime clustering occurs by caustics. 



number, but an additional dimensionless parameter, re, 
defined below. We infer that strong clustering can occur 
at large Stokes numbers. Two distinct mechanisms com- 
pete (clustering onto fractal sets versus clustering onto 
caustics in an otherwise homogeneous background) and 
dominate in different regions of the parameter space. 

We model the particles suspended in the fluid flow by 
the equation of motion 



l{u{r,t) - r) 



(3) 



where r = (ri,r2,r^) denotes the position of a particle. 
Eq. @ is appropriate for non-interacting spherical par- 
ticles when the Reynolds number of the flow referred to 
the particle diameter is small. It is assumed that the ra- 
dius of the particle and the molecular mean free path of 
the fluid are sufficiently small. Stokes's formula gives the 
damping rate 7 = Qiraptv / m where v, pt are respectively 
the kinematic viscosity and density of the fluid, and a, 
m are the radius and mass of the particle. Effects due 
to the inertia of the displaced fluid are neglected. This 
is justified when the density of the suspended particles is 
large compared to that of the fluid. We also assume that 
Brownian diffusion of the particles is negligible. 

We now discuss the dimensionless parameters of the 
problem: the velocity field is assumed to be characterised 
by its typical velocity u = (it 2 ), by a correlation length 
£ and a correlation time r. In addition, the interaction of 
the fluid with the particles is determined by the damp- 
ing rate 7. From these four quantities we can form two 
independent dimensionless groups: a dimensionless ve- 
locity re = ur/t;, and the dimensionless damping u> = jt 
(so that St = u >~ j . The parameter re has been termed 
'Kubo number' |L7}. It has not been considered before in 
this context. We argue that it cannot be large if u(r,t) 
is to be a satisfactory model for a solution of the Navier- 
Stokes equations: t < £/u since disturbances in the fluid 
velocity field u(r,t) are transported by u(r,t) itself. 

Consider now the particular case of fully-developed 
turbulence. In this case, the velocity field exhibits a 
power-law energy spectrum, with upper and lower cut- 
offs 0]. The smaller length scale is the Kolmogorov 
length, which is the size of the smallest vortices gener- 
ated by the turbulence. It is given by (^ 3 /e) 1 / 4 , where 
e is the rate of dissipation per unit mass of fluid. The 
Kolmogorov length corresponds to the correlation length 
£ in our theory. The corresponding typical velocity u and 
correlation time r are also determined solely by the same 
two parameters, e and v, implying that k ~ 1 for fully 
developed turbulence. In other situations re can be small. 

We now turn to a summary of our results and outline 
how they were derived (details will be published else- 
where). Linearising the equations of motion JSJ gives 

dp = —jSp+F(t)5r, 5r = Sp/m (4) 

where p = mr is the particle momentum and F(t) is 
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matrix of force gradients: 



7m^(r(t),t) 



(5) 



We take three trajectories displaced relative to a ref- 
erence trajectory by small increments (Sr^,Sp^), with 
fi = 1,2,3. We introduce a triplet of orthogonal unit 
vectors n„(t) such that ni(i) is oriented along 6ri(t), 
and n2(i) lies in the plane spanned by (Sri(t), Sr2(tj). 
This determines ^(t) up to a sign which is fixed by 
requiring continuity. We write n„(i) = O(i)n„(0) and 
Sp^t) = R(i) 5r^(t) where O is an orthogonal and R a 
general 3x3 matrix. We define the elements of F and R 
transformed to the moving basis: 

F^(t) = n^(t)-F(t)n„(t), ^(t)=n M (t)-R(t)n v (t) (6) 

and find the following equation of motion for R' 

R = - 7 R - — R' 2 + [R, O+O] + F' . (7) 

TO 

The elements of + are given by 
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We find that the Lyapunov exponents are equal to the 
long-time average of the diagonal elements of R' 

A : = (i?' n )/TO, A 2 = (R' 22 )/m, A 3 = (R' 33 )/m. (9) 

Eqs. Q and JHJ for R' can be simplified when the cor- 
relation time of the velocity field is sufficiently short, 
w<l, assuming that the amplitude of the random force 
is sufficiently small, n <C 1. In this limit F' behaves 
as a white-noise signal, and 10 reduces to a system of 
Langevin equations. We label the dynamical variables 
by a single index i = 3(/x — 1) + v and scale the Langevin 
equations for B! { to dimensionless form 

9 9 

dxi = - (xt + e ^2 ^2 1 ./< •' .••''< + dw » ( 10 ) 

j=l k=l 



Here t' — -ft, Xi = y/j/DiR^, and (dwidwj) — 2Dijdt'. 
The elements Dij of the diffusion matrix D are given by 



three Lyapunov exponents from the expectation values 
of variables in a system of Langevin equations. Earlier 
work has obtained the largest Lyapunov exponents for 
various problems using Langevin equations [l3 |2fJ, |2lJ . 

The elements of D are determined by the fluctuations 
of the velocity field. We assume that the latter is in- 
compressible, but for reasons explained below we add a 
small compressible component: u = VA A + VSAq. The 
fields A^(r,t), fj, = 1,2,3 are taken to be homogeneous 
in space and time, and isotropic in space. Their correla- 
tions are determined by {A^{r + R,t + t)A u (ro, to)) = 
<W C(|r - r'|/£, \t - t'\/r). The field SA a is statistically 
independent of A^, has the same correlation function, 
and in the end the limit SAq — ► is taken. 

The Langevin equations IjlQI) are equivalent to a 
Fokker-Planck equation whose stationary solution P(x) 
determines the Lyapunov exponents. In the limit of 
e — ► the latter is Gaussian 



Po(x) ex exp(— |cc • D 1 x) = exp[— ( E ) o( a; )] 



(13) 



This suggests transforming the Fokker-Planck opera- 
tor so that its e = limit is transformed into a 
harmonic oscillator. This is achieved by introducing 
Q(x) = exp[$ (x)/2]P(x). The steady-state Fokker- 
Planck equation can be written as (Ho + eHi)\Q) = 0, 
where we have represented the function Q(x) by a 'ket 
vector' \Q). The operator Ho is the Hamiltonian for nine 
uncoupled harmonic oscillators 



H = - ^ "'t&i 

i=l 



(14) 



where the and a,i are, respectively, the creation and 
annihilation operators for the degree of freedom labelled 
by i (satisfying [dj,a^] = SijI). The non-Hcrmitean per- 
turbation Hi can be expressed in terms of the eigenvalues 
LOi of D and the elements of an orthogonal matrix J 
satisfying D = Jf2J _1 , with ft = diag(wi): 

#i = ^Bglaf (a} + + a k ) 

ijk 



- ^UjLUk/uJi y^VmnJilJmjJnk ■ (15) 
Imn 



Di3=\ 



dt(Fl(t)F'(0)) 



(11) 



The coefficients V^ k are determined by the 2nd and 3rd 
terms on the rhs of JJJ. The dimensionless parameter 



e = D{{ 2 /( mi 3 / 2 ) ~ ku- 1 ' 2 



(12) 



is a measure of the inertia of the particles: it is propor- 
tional to a and therefore to to 1 / 3 . Thus we obtain all 



Regularisation is needed since one eigenvalue vanishes in 
the limit of SAq — > 0. We determine \Q) by perturba- 
tion theory in e. Given \Q), the Lyapunov exponents are 
obtained as Ai = je(xi), X 2 = 7e(xs), A3 = je(xd), and 



(Xi) 



1 



(*o|Q) 



J2JvV^(^o\a 3 +a+\Q) (16) 



where |$o) denotes the null eigenvector of Ho- From l|ltj|) 
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we obtain series expansions in the form 

A1/7 = 3e 2 - 29e 4 + 564e 6 

-14977e 8 + 488784e 10 - 18670570e 12 + ■ ■ ■ 
A 2 / 7 = 8e 4 - 459/2 e 6 + 14281/2 e 8 (17) 

-757273/3 e 10 + 361653709/36 e 12 + • • • 
A3/7 = -3e 2 -9e 4 - 789/2 e 6 - 5787/2 e 8 

-895169/3 e 10 - 101637719/36 e 12 + • • • . 

Note that only even powers of e contribute, and that 
all coefficients are rational numbers. Eq. (|17|l is the 
main result of this letter. The expansion is valid in the 
underdamped limit lu <C 1 when k«1. 

The coefficients in (I17|l exhibit rapid growth typical of 
an asymptotic series |22| . We have attempted to sum the 
series ifTTjl using Pade-Borel summation 22] 

, W* Jj) 

A,-/ 7 ~Re / dte-^^ 2 ' (18) 
Jc 1=1 

where Cj are the coefficients of <(T7|) and Z max = 7 is 
the number of nonzero coefficients available for each \j. 
The sum in the integrand is approximated by Pade ap- 
proximants |2^| of order n, namely P™ or P™ +1 with 
n < [imax/2]. The integration path in H18[) is taken to be 
a ray in the upper right quadrant in the complex plane. 

Results of Pade-Borel summations of the series for Xj 
are shown in Fig. ^1 and converge to results of nu- 
merical simulations provided e is not too large. For A2 
numerical evidence indicates the presence of additional 
non-analytical contributions not captured by the Pade- 
Borel summation. The results of Fig. la allow us to 
determine the quantity A defined in eq. @. The re- 
sult is shown in Fig. lb. We find that A is maximal 
for e « 0.21 and positive (indicating clustering onto a 
fractal set) for < e < 0.33. The red line in Fig. QJi, 
e ~ kuj^ 1 / 2 = const., indicates schematically where A is 
zero. Above the red line A is always positive, but tends 
to zero for small e as A = 10e 2 ~ k 2 /lu. In the limit 
of e — > the dynamics becomes advective (despite being 
underdamped): to lowest order in e our results coincide 
with those for purely advective flow ||. 

We now turn to the rate J of caustic formation. It is 
the rate at which 5V(t) = (Sri(t) A <5r 2 (i)) ■ 5r 3 (t) goes 
through zero. Since 5p^ typically remain bounded, caus- 
tics correspond to instances where the elements of the 
third column of R' go to —00 and reappear at 00. The 
rate at which these events occur is given by the escape 
rate of the Langevin process l|10|) to infinity. It is ex- 
pected 0] to have a non-analytic dependence on e, of 
the form C exp(—S/e 2 ), as demonstrated in Fig. lc. In 
this panel, J/7 is compared to — (Ai + A 2 + A 3 )/7. We 
see that caustics are very rare when e -C 1, but frequent 
when e is large and they are the only clustering mecha- 
nism when e > 0.33. 



Finally, we comment on the relation between our re- 
sults and earlier works (cited above), which suggest that 
clustering only occurs for w k 1 (with the value of k un- 
specified). It must be emphasised that the earlier quan- 
titative theoretical results on clustering are confined to 
the overdamped limit u> 3> 1, where inertial effects are 
small: for purely advective flow there is no clustering 
(A = and J = 0). Inertial effects were incorporated 
by Elperin and others 0, Q, 0] , using Maxey's pertur- 
bative correction to the velocity. Their results are valid 
only for the limit lu 3> 1, and are distinct from our series 
expansions (|17|l : this is most easily seen by calculating 
corrections to A in this overdamped limit. We find that 
A ~ k 2 /uj 2 implying that clustering effects are small in 
this regime. In the underdamped regime, by contrast, we 
obtained A ~ k 2 /lu which can be of order unity. 

The results of this letter are summarised schematically 
in figure Id. First, at small k, strong clustering occurs in 
the region indicated, above the line e ~ klu^ 1 / 2 = 0.33. 
Second, since the dimension deficit A is positive in this 
regime, the reasoning of Sommerer and Ott |14| indicates 
that the particles cluster on a fractal. Third, as e — * 
the dynamics becomes advective. In this limit the dimen- 
sion deficit A and the rate of caustic formation J vanish: 
particles advected in an incompressible flow remain uni- 
formly distributed. Fourth, when e > 0.33 we find that 
the dimension deficit A is negative implying that do not 
lie on a fractal. They are however not homogenously dis- 
tributed: in this regime particles cluster because they are 
brought into close contact by caustics. 
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